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ABSTRACT 


The stream function, vorticity form of momentum equations together 
with the equation for conservation of energy are solved numerically using the 
Spectral Element Method for two and three dimensional natural convective 
flows. For the two dimensional situation, two physical problems, such ais, 
flow and heat transfer in a buoyancy driven cavity and natural convection 
in a closed cavity with discrete heat sources have been considered. For the 
three dimensional situation, computations have been performed for natural 
convection heat transfer from rectangular vertical fin arrays. The streamlines, 
the vorticity distribution, and the temperature fields for various values of 
Rayleigh and Prandtl numbers are obtained. Subsequently, the influence 
of different geometrical parameters, Rayleigh and Prandtl numbers on the 
flow and heat transfer is investigated. An attempt has been undertaken to 
compare the predicted results with the results reported in open literature by 
other investigators. 

The above said problems find applications in cooling of electronic equip- 
ments and integrated circuit packages. 
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NOMENCLATURE 


a,b 

Ar 

[C*], [C] 
{e*}, {e} 

f 

{f} 

g 

Gr 

H 

I 

L 

Ny 

Nu 

P 

P 

Pr 

Ra 

t 

T 

u,v 

U,V 

x,y,z 

X,Y,Z 


Elemental location 
Aspect ratio 

Elemental left side coefficient matrix 

Elemental right side coefficient matrix 

Number of elements in x direction 

Number of elements in y direction 

Forcing function in Helmholtz equation 

Array of forcing function at the grid points 

Acceleration due to gravity 

Grashof number 

Height of the cavity 

functional 

Length of cavity 

Number of Gauss-Lobatto Chebyshev points in x direction 
Number of Gauss-Lobatto Chebyshev points in y direction 
Nusselt number 
Pressure 

Nondimensional pressure 
Prandtl number 
Rayleigh number 
time 

Temperature 
Velocity components 
Nondimensional velocity components 
Cartesian coordinates 
Nondimensional cartesian coordinates 




Nomenclature 


XU 


Greek symbols 

a Thermal diffusivity 

j3 Cofficient of volumetric expasion 

r Diffusion coefficient 

0 Non dimensional temperature 

A Dummy parameter 

A A constant in general transport equation 
p Density 

V Kinematic viscosity 

r Non dimensional time 

At Non dimensional time step 

(f) Unknown function in general transport equation 
rj) Stream function 

(jj Vorticity 

V Nabla operator 
Laplacian operator 

Subscripts 

c cold 

h hot 

1 matrix index 

j,k,p,q matrix indices 

00 Reference state 

Superscripts 

1 elemental parameter 
n time step 



Chapter 1 


Introduction 


1.1 A Review of the Numerical Methods 

The techniques that are most commonly used to solve flow problems in computational 
fluid dynamics are the finite difference method, the finite element method, and the 
spectral method. 

The finite difference method has been around the longest and is still popular 
because of its simplicity. 

The finite difference method generally provides pointwise approximations to par- 
tial differential equations. Derivatives in the governing partial differential equations 
are replaced by equivalent finite difference expressions which involve the values of 
the solution at discrete grid points of the domain. Finite-difference method is easy 
to implement but it suffers from several disadvantages like low finite-order accuracy 
in evaluation of the derivatives, the difficulty in handling irregular geometries, the 
difficulty in imposing unusual boundary conditions along arbitrary boundaries, and 
the inability to employ nonuniform grids with ease. 

Finite element method has been widely applied to many applications related to 
mechanical engineering, especially in solid mechanics, because of its extreme geomet- 
rical flexibility. These methods are also used to solve heat transfer and fluid flow 
problems. 

Finite element methods provide piecewise, or regional approximation to partial 
differential equations. Various methods of this class exist; all requiring that an inte- 
gral representation of a partial differential equation be constructed. Classical finite 
element methods for structural mechanics are based on variational principles. For 
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transient diffusion and conduction and for convective heat transfer processes more 
general approaches, such as a method of weighted residuals, are used. Finite element 
method is a low order h-type weighted residual method. 

The principle advantage of the finite element approach for convective heat transfer 
is in the ability to handle irregular geometries. It generally requires fewer node points 
for a given accuracy as compared with finite difference methods. The drawback of 
the finite element method is that it is a low order method. Like finite difference 
methods, the finite element method converges slowly as the number of the elements 
is increased. 

Finite difference and finite element methods are finite order accurate, i.e., the 
error of approximation behaves like (1/iV)’”, where N is the number of nodal points 
and m the order of the scheme. 

The spectral method is a high order weighted residual method [l].It is also known 
as the p-type weighted residual method.Unlike the finite element method, the trial 
functions of the spectral method are global, so it is usually applied to regular geometry 
problems with smooth solutions. The spectral method is known for achieving high 
accuracy. For equivalent number of grid points, spectral methods provide excellent 
accuracy compared with either finite-difference or finite-element methods. It has 
been estimated that in order to obtain 1% accuracy for a wave-like solution, 40 grid 
points per wavelength are needed for the second-order finite difference method,but 
only 3.5 collocation points per wavelength are needed for the Chebyshev spectral 
method. Therefore the spectral method is about 10 times more efficient than the 
finite difference method in storage aspect alone ( in 1-D ).Peyret and Taylor [10] 
compared the finite difference method with the spectral method and found that the 
spectral method was at least 10 times faster than the finite difference method for 
the same accuracy in the vortex problem they solved. Spectral methods display 
exponential convergence rates forincreasingN. The convergence rate associated with 
the spectral methods depend only on the smoothness of the approximated function 
within the computational domain and the approximation error decreases faster than 
any finite power of 1 /N, where N is the number of nodal points in the domain. 

The basic idea of spectral method is to express the dependent variables in the 
form of a truncated expansion of orthogonal eigen-functions of the Sturm-Liouville 
problem. The coefficients of this expansion are determined using either a Galerkin 

(MJL. 

or a Collocation formulation. The trial functions which^commonly used in spectral 
expansion include the Chebyshev, Legendre, and Laguerre polynomials for problems 
involving nonperiodic boundary conditions. For periodic problems the global inter- 
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polant is a Fourier series. 

However, the spectral method has its own limitations. For instance, the grid 
interval is not uniform and it decreases drastically toward the boundaries of the com- 
putational domain. This imposes a severe limitation on the maximum time step 
when solving initial value problems via an explicit time marching scheme. Therefore, 
we often use implicit methods in conjuction with spectral algorithms. Unfortunately, 
these implicit schemes are not very efficient either, because the matrices in the spectral 
method are full matrices, unlike the narrow banded matrices of the finite difference 
and even the finite element methods. The main weakness of the spectral method 
relativite to FDM and FEM is its inflexibility to be adaptable to irregular com- 
putational domain. The solution of the physical problems should be smooth and the 
solution domain should be rectangular-like; otherwise spectral accuracy cannot be 
obtained. The appearance of ‘spectral element method’ that combines the generality 
of FEM and the accuracy of SM, has overcome the problem of restriction on the 
geometry. 


1.2 Spectral-Element Method 

The spectral element method (SEM) is a synthesis of the spectral and the finite 
element method [2,3]. The spectral element method is also called the ‘p-type finite 
element’ method, or the h-p type weighted residual method. Like the spectral method, 
it uses high order polynomials as trial functions,but like the finite element method, 
it decomposes the computational domain into many elements and defines local trial 
functions. The hybrid character of the spectral element method enables it to overcome 
the shortcomings of both the spectral method and the finite element method but still 
retain their advantages. Since the trial functions of the spectral element method 
are local, it can handle complex geometry easily. On the other hand, it is still a 
high order weighted residual method, so the exponential convergence rate is achieved 
as the degree of the polynomials in each element is increased. However, eis with 
the conventional spectral method, the spectral accuracy can only be achieved when 
the solution possesses no singularities. In the spectral element approach, the spectral 
approximation is applied on each of the subdomains. Therefore, for those subdomains 
where the solutions are not smooth, the high order accuracy of spectral methods will be 
destroyed. The main difference between the spectral element method and the spectral 
multi-domain method is that the (7° and boundary conditions at the interface of 
the elements have to be explicitly enforced by the spectral multi-domain method. The 
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spectral element method, by contrast, uses the variational principle to guarantee 
and (weakly) continuity at the interface, which results in a much simpler and more 
natural approach than the non-variational method; therefore, parallel algorithms can 
be conventiently implemented [3]. 


Chapter 2 


Mathematical Formulation 


2,1 Introduction 

The spectral element method can be applied to simple (e.g., rectangular, square do- 
mains) as well as arbitrary (curved) geometries. However, for complicated geometries 
the technique can be applied only after coordinate transformation. The complicated 
geometry is first changed into a simple geometry (rectangular or square) by some 
isoparametric transformation and the technique is then applied. Any complex ge- 
omet^ in two dimensional space can be broken into multiple elements of standard 
rectangular or quadrilateral shapes. 

In the spectral element discretization, the computational domain is broken into a 
number of elements, and within each element the dependent variable is represented 
as a high-order Lagrangian interpolant, in terms of Chebyshev polynomials, the coef- 
ficients of which are related to the function values at the Gauss-Lobatto Chebyshev 
collocation points. 

The governing equations are partially discretized using finite difference approxi- 
mation for the time derivative. The convective and source terms are treated explicitly 
and the diffusion term (elliptic contribution) imphcitly. 
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Table 2.1: Expressions for (f), and in general transport equation (dimensional 
form) 


Equation 

4> 

F 

Sff, 

Continuity 

1 

0 

0 

X momentum 

U 


-|^ + {Poo - p)9x 

y momentum 

V 


“ ^ d" {poo ~ p)9y 

Energy 

T 

Jc_ 

Cp 



Table 2.2: Expressions for and in general transport equation (nondimensional 
form) 


Equation 

4> 

F 



Continuity 

1 

0 


0 

X momentum 

U 

1 

HI 

dx 

+ Su 

y momentum 

V 

1 

H| 

By 

+ Sv 

Energy 

d 

1 

Pr 


Sx 


2.2 Governing Equations for Laminar Incompress- 
ible Flow 


The general transport equation for the flow property <j> in cartesian coordinates can 
be written as 


d4> fTT^^ 




dX^ dY^ 


+ Sa 


(2.1) 


Where is a source term , F i« a diffusion coefficient. The parameter <f> represents 
l,u,v,T,w,V’,or^. A summary of the and F^ applicable to each (f> is given 
in Table 2.1. 

Writing the general transport equation in a dimensionless form facilitates general- 
isation to embody a class of problems. In the usual dimensionless form of the general 
transport equation the parameters are defined in a way as given in Table 2.2. 
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Table 2.3: Expressions for <f>, p, and S 4 , {il> — u>)approach 


Function 

<t> 

p 

r 

S4, 

Temperature 

e 

1 

1 

Pr 

St 

Vorticity 

w 

1 

1 


Stream function 

V’ 

0 

1 

to 


If the governing dimensionless equations for the conservation of mass , momentum 
and energy are formulated in terms of vorticity , stream function and temperature then 
the parameters used in general transport equation (2.1) are defined in Table 2.3. 


2.3 Time Discretization 


Finite-difference approximaton is used to handle time dependent term. Both explicit 
and implicit methods have been used to a variety of flow situations. Generally, the 
imphcit methods have been recommended because of their favourable stability prop- 
erties. The explicit schemes for parabolic unsteady partial differential equations axe 
required to satisfy certain stability requirements, which impose a restriction on the 
time step. Since the implicit methods are generally free from these restrictions, it may 
be possible to use a relatively larger time step than what is permissible in explicit 
schemes. However, an implicit method has its own limitations. An implicit scheme 
requires a relatively large number of arithmatic operations than an explicit scheme. 
For a coupled non-linear partial differential equation set, the impHdt formulation will 
be even more time consuming because an iterative scheme is required to take care of 
the non-liearities and the coupling in the governing equations. 

Patera [2,3] uses an Adams-Bashforth Crank-Nicolson (ABCN) scheme with a 
spectral element spatial discretization. 

To solve the general transport equation [Eqn (2.1)], a semi-implicit, time-marching 
scheme is employed for time discretization, wherein all the linear terms (continuity, 
pressure and diffusion terms) are treated impficitly and the non-linear convective 
terms explicitly. The time-discretized form of equation (2.1) is given by 



p [dX^ dY^\ ^ P 


( 2 . 2 ) 
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which can be rearranged in the form 

\M^) "^V^j ~TK^ 


TArf p 


1 

T/P 




n cn 

- -f (2.3) 


or 


«>.. + «S„-AV=/ (2.4) 

This is the well known Helmholtz (Second-order elliptic) equation. The function <f> 
is unknown function to be computed and ‘ f ’ is the known forcing function. Here n 
refer to the known values at time r, while all other values, n + 1, are unknown values 
at time r + At. The parameter can take different positive values. In the limiting 
case when — > 0, equaton (2.4) reduces to Poisson equation, which too can be 
solved using the Helmholtz solver. 


2.4 Spectral Element Spatial Discretization 

In order to solve the Helmholtz equation [Eqn. (2.4)] by spectral element method we 
use its equavalent variational formulation (as in the finite-element technique) which 
amounts to the maximization of the functional 

/ = // + (4',?} - ^4? - 4f dxdy (2.5) 

The integration in equation (2.5) spans over the entire computational domain. 

When the domain is broken up into M elements, the functional will simply be the 
summation of the contributions from the individual elements, namely 

M 

/ = EC (2.6) 

1 = 1 

For the element the functional, is given by 

r=Jj [- 5 {(«i.'y + (V)“} - - 4'r dxdy (2.7) 

The intregration in equation (2.7) spans over the area of the element , defined 
as the rectangle between [aj,,6j.] and (In this study we use only rectangular 

elements.). 
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Within the element the variable <f>{x,y) is represented cis a Lagrangian inter- 
polant, in terms of Chebyshev polynomials, through + 1) Gauss-Lobatto 

Chebyshev points given by 


x} = 

= cos(— ) 

^ X 

j — 0) 1) 2, 


(2.8) 

yi = 

vk 

- cosij^) 

A: = 0,1,2, 

TV* 

(2.9) 


y 


The Lagrangian interpolant of a variable (f>{x, y) in the element (of dimensions 
LI X X‘ ) is given by 


K 

j=0 k=0 


( 2 . 10 ) 


where are the values of <l> at the point {xj,yk) and {x,y) is the local coordinate 
system normalised so that the rectangular elements lie between x = ±l,y = ±1. 
The local element coordinate system is therefore defined as 

x‘ = -^(x-4)-l ( 2 . 11 ) 


y' = -^(y - O - 1 (2-12) 

^y 

where L\.{= — o^) and Lj(= hy — are the lengths of the element in the x and y 

directions respectively. 

A local coordinate system (i.e., a coordinate system in the element) proves to be 
convenient in the use of the interpolation functions.The global coordinate x (which is 
used to describe the governing differential equations of the problem) can be retrieved 
from the local coordinate by inverting the linear transformation given by equations 
(2.11) and (2.12). 

The X direction interpolation functions hm{x) are order local Lagrangian 
interpolants in terms of Chebyshev polynomials satisfying 


( 213 ) 

within the element (at the local coordinates); 6mn is the Kronecker-delta symbol and 
x„ denotes the x coordinate of the Gauss-Lobatto point in the x direction, ijj, 
are identically zero outside the element. 
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The interpolation functions hjj,(x^)can then be expressed as 



n=0 

(2.14) 

where the T„ 
als Tn{x) are 

are the order Chebyshev polynomials.The Chebyshev polynomi- 


Tn{cosd) = cos{n9) 

(2.15) 

or 

Tn{x) = cos{ncos~^x) 

(2.16) 

and therefore 

Tn{Xj) = COs{ ) 

(2.17) 


also in equation (2.14), the integer functions C are 

Ck = 1, for or iV’ 

= 2, for A: = 0 or iV’ 

We approximate the variable <j) in the element by 

<f>^ 4 >n <l>jkhj{x)hk{y) 

j=Q k=0 


(2.18) 

(2.19) 


( 2 . 20 ) 


where hk{y) are interpolation functions similar to hm{x) explained above, except that 
it is for Ny points. 

We can compute the maximum error in the domain between the exact solution 
and the calculated solution (j>]^ using the norm 


11= max\<j){x,y) - 


( 2 . 21 ) 


for 


-l<x<l -l<y<l 
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This measure of difference is called supmetric.Note that the supmetricis a real number 
indicative of error. This error reduces very fast as the number of points is increased 
in a spectral method. 

With these particular collocation points, the Lagrangian interpolant for satisfies 

IK NiNl oo (2.22) 

for all finite k. This property is called ‘exponential’ or ‘infinite-order’ accuracy of 
interpolation. 

Other choices of collocation points have the property (2.22), however, (2.8) and 
(2.9) seems a reasonable choice given the good approximation properties of Chebyshev 
polynomials as well as the existence of a fast transform for Chebyshev series. 


2.5 Formation of Elemental Equations 

The forcing function f, in Helmholtz equation (2.4) is also interpolated in the same 
fashion as is interpolated.For element 


noting that 


Ni Ni 
j=0 Jt=0 


2 1 

n=0 


(2.23) 


(2.24) 


Therefore 


JV* 

2 1 




N* N' N* N' 

T„(a)T„(ir) 


and 


Ni Ni Ni Ni 


5’) = E E E E 

j=z0k=0n-0m=0 

Tm{yi)Tmm 


(2.25) 


(2.26) 


(2.27) 
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In order to apply spectral element method involving Chebyshev polynomials each 
element should be mapped (using linear transformation equations (2.11) and (2.12)) 
onto a domain spanning from -1 to -|-1 in x and y directions respectively, because 
Chebyshev polynomials itself are defined over this span. Whenever a differentiation 
or integration is performed in the transformed co-ordinate (x, y), it should be properly 
scaled. 

Since 


X = 


— (x -a^)-l 
■L/x 


We have 



( 2 . 28 ) 


( 2 . 29 ) 


The differentials with respect to x , <f>x and <f)s in the two coordinate systems are 
related by 


d<j) _ d(f>dx _ 2 d(j) 
dx dx dx Lx dx 


( 2 . 30 ) 


d(j) _d<j)dy _ 2 34 
dy dy dy Ly dy 


( 2 . 31 ) 


where and ^ are scale factors. The derivatives || and %, can be obtained as 


dy 


Ni Ni 

'/A*.! __ 

j=0k=0 

After substitution this equation reduces to 


34^ ^ dk'{x') ■ , 


34* 

3x 


Ni Ni 


Nl 






Nl 


j=0 k=0 


jk j^i n,n 


■UA) 


tE 


d(x<) 


The resulting expression may be written as 


9 zt ■^•1 . 1 

= IT MK § S n?o 5, CfiAC, 




( 2 . 32 ) 


( 2 . 33 ) 


( 2 . 34 ) 
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Similarly 




at; n{ w; 


77 E E E E C.C 

j=0 Ar=0 n=0 mnO 




(2.35) 


We now insert the interpolants (2.26), (2.27), (2.34), and (2.35) into the func- 
tional (2.7), perform the resulting integrals, and require stationarity with respect to 
variations in the i.e., require that the variation of with respect to the nodal 
values vanish : 


0 

II 

(2.36) 

Remember that 



(2.37) 

Then the elemental algebraic system of equations will simply be 


^jklm^lm ^jklmflm ^Im 

(2.38) 

^jklm = ^jklm ~ 

(2.39) 

Here 


r t 1 r * 1 

^jkim Tj, (iV‘iV‘)2^ LiiNiNif 

(2.40) 

r t r t -1 

pt „ -^x^y J- B^^hhy 

^jklm ^ (N'N^y 

(2.41) 



4 


Ni Ni 

EE 


n=0 m=0 



Tn{x'j)Tmix])anm 


(2.42) 


A I 

p*,» — V' _ _ 

CiC,hii.cA 




(2.43) 



Mathematical Formulation 


u 


and 


Where 


Also 


^nm 


L 


0 , 

nm 


-1 dx' dx' 


dx' 


[J\(n—m)l 2 \ «^|(n+m)/ 2 |]) 


nFm odd 
n + m even 


(2.44) 




1 

2p-l’ 


h > 1 
k = 0 


(2.45) 


j ^ TnTmdx' 

0 , 

1 1 

1 — (n + m)2 1 — (n — m)^ ’ 


n + m odd 
n-\-m even 


(2.46) 


2.6 Formation of Global Matrix Equations 

Once the system of discrete equations for each element is obtained, with proper 
boundary conditions imposed for the elements having one or more sides lying on 
the boundary of the computational domain, the system of global matrices is then 
constructed by ‘direct stiffness’ summation, adding at the boundary nodes the con- 
tributions from the corresponding elements. 

Owing to the difficulty in direct implementation with four dimensional arrays [B] 
and [C] in equation (2.38), this elemental matrix equation is transformed into a more 
conventional form by combining the two indices j and k to give a third index p, and 
the indices 1 and m to give a third index q by the following equations 

p={N, + l).j + k (2.47) 


where 

i = 0,l,...,iVa„ fc = 0,l,...,iVj, 


SO 


p = 0, 1, ..., {Nx + -f 1) — 1 

(2.48) 

and 


q = (iVy + 1).Z + m 

(2,49) 
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( 0 , 5 ) 

( 1 . 5 ) 

( 2 . 5 ) 

( 3 . 5 ) 

( 4 . 5 ) 


Figure 2.1; Original numbering of grid points in an element 

/ = m = 0,l,...,JVj, 

or 

5 = 0,1,...,(JV. + 1).(W, + 1)-1 (2.50) 

when this is done, equation (2.38) is transformed to convetional matrix form 

= B',/‘ = e' (2.51) 

Any grid point in an element, which was specified by two indices in equation (2.38), 
is now specified by a single index, say p following the above transformation. This 
changed numbering system for Nx = 4 and Ny = 5 is illustrated in Figure 2.1. The x 
direction is taken to be vertically downward and y direction horizontally rightward. 

Based on the problem requirements, the computational domain may consist of a 
number of elements. The following remarks can be made on the comparison of the 
multi-elements and single-element case, keeping the total number of grid points in 
the computational domain exactly the same in both the cases [1 1] . 

• The single-element method will give a higher order approximation. 

• The multi-element method will give greater geometric flexibility. 

• The global matrix for the single-element case will be full while that of the 
multi-elements case will be relatively sparse. If this sparsity can be exploited 
by the solution procedure the multi-element method can be made more efficient 
in terms of computer time and memory. 
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Figure 2.2: Changed numbering of grid points in an element 

The higher order approximation in single element case results in higher accuracy 
compared to multiple elements. The accuracy depends upon the order of the La- 
grangian interpolants and not on the total number of grid points. In multi-element 
case, to keep the total number of grid points same in the domain, the number of 
points per element get reduced and consequently the function will be approximated 
by a lower order interpolant, thereby causing deteriotion in the solution. 

The conclusion drawn from the above discussion is that the main advantage of 
using multiple elements over single element lies in the greater geometric flexibility 
afforded by the former. 

The total number of grid points in the domain in case of single element will be 
{Nx + l)iNy + 1). The corresponding coefficient matrix [C] will be of order {Nx + 
l)(Ny-\-l) X (Nx + l){Ny + 1). If a rectangular domain consists of elements in the 
X direction and Cy elements in the y direction, then the total number of grid points in 
the domain for Cx x Cy elements will be l)(iVy.ej,-}-l) and the global coefficient 

matrix will be of the order {Nx-tx + l)(iVj,.ej, + 1) x {Nx-Cx + l)(iVj,.ej, + 1). 

The global grid points numbering for a system of 2 elements in x direction and 2 
elements in y direction, i.e., = 2 and ey = 2 and Nx = 4:, Ny = 5 is shown in Figure 

2.3. Listed in parenthesis are the element numbers, e.g., (k,l) represents the 
element in the domain. 

The global numbering of grid points of {k, element is as follows: 

• gird point numbering on left face of {k, ly^ element is given by 
{Ny.ex + l){(fc - l)Nx + *} + [Ny.l +1)-Ny 
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• grid point numbering on right face of {k, /)** element is given by 
{Ny.e, + 1){(A: -l)N, + i} + {Ny.l + 1) 


Every grid point in the computational domain ha^ two numbers associated with 
it: a local and a global number, as shown in Figures 2.2 and 2.3. The local number 
corresponds to its position in its element and the global number corresponds to its 
position in the whole computational domain. 

Construction of global matrix, from elemental matrices, [C], can be ac- 

complished as follows: 

For every point that lies in the interior of an element, the coefficients of the el- 
emental matrix [C] in a locally numbered row and colunm are transferred to the 
corresponding globally numbered row and column in the global matrix with- 

out any modification in nodal values. Same procedure is applied for elemental matrix 
{ e } and global matrix In doing so, a proper correspondence is made 

between the local and global number of the grid point in the domain. 

However, at the intersection of the elemental boundaries in the domain, the equa- 
tions corresponding to grid points at the interface of the two elements are obtained 
by adding the two equations corresponding to these grid points in the two elemental 
sets, 

The global equations then are 

( 2 - 52 ) 

where the vector incorporates all the unknown (f>im of all the elements and the 

indices pg and qg represent the global indices, given by 

Pg = 1, ..., (iV^.ex -t- l).(iVj,.ej, -b 1) (2.53) 

qg — 1, •••5 (^Nx-Cx “b l)-(7Vy.ey "b 1) 

The inversion of the global system matrix is carried out directly ( LU decompo- 
sition ) in a preprocessing stage at the beginning of the simulation. Thereafter, only 
matrix multiplications are performed at each time step,which is very convenient re- 
garding processing time. 
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2.7 Imposition of Boundary Conditions 

2.7.1 Imposition of Dirichlet Boundary Conditions 

Dirichlet boundary conditions are imposed by matrix condensation. The rows and 
columns corresponding to boundary points are eliminated from the system matrix. 

There is another way of imposing these boundary conditions. The diagonal ele- 
ments of in equation(2.52) corresponding to the boundary points are equated 

to a very large number and the corresponding elements in the right side of are 

set equal to the specified boundary value multiphed by the same large number. Then 
the resulting matrix equation is solved to obtain the same solution. This method is 
more expensive as it solves the full matrix equation with all points rather than just 
the points where the values are unknown. However, it has the advantage of simplicity 
of implementation. 

2.7.2 Imposition of Neumann Boundary Conditions 

Mathematical treatment of the spectral element technique is such that the zero Neu- 
mann conditions are naturally imposed and no modification of the matrix equations 
(2.38) is needed to impose these boundary conditions. However, non-zero Neumann 
conditions are imposed by modifying the functional (2.5). 

To illustrate the procedure to impose the non-zero Neumann conditions, let us 
assume that 

|^ = A(j() at 1 = 4 (2.54) 

and the boundary of this element coincides with the boundary of the domain. 

According to variational formulation, solution of the differential equation [Eq.(2.4)] 
subject to given boundary condition [Eq.(2.54)] is equivalent to maximization of the 
functional 

m = // f- + (4)'} - H 

where he is the boundary term corresponding to the non-zero Neumann condition 
[Eq.(2.54)]. For each element i which has that boundary condition in the domain, 
Eq.(2.7) can be written as 

m = /J' + (4?) - - «■] d^dy + (2.56) 
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I I . {-My)4>}dy\x=o 

Jal Jal 


i.e., 


r - r -I- p 

intertor * -^ax 


(2.57) 


Because of the additional term in equation (2.57) compared to equation (2.5), the 
final elemental matrix equation [Eq.(2.38)] corresponding to the grid points on this 
boundary will have some additional terms in the right hand side. These terms can 
be determined as follows: 

The interpolant of variable (f> on this boundary of element is a function of y 
only, i.e.. 






similarly 


k=0 n=0 


jyy 2 1 

Tlfo) = E 7r7rT„{y))T^(v‘) 

i=0 mzzO 


(2.58) 


Therefore 


II 


(2.59) 


(2.60) 


- P{A<l>%^^idy 

- j_^A{y).4>{y)}^dy 

-f /.Ea4 E 


N', 


E 7?7rTM)Uf) 

n=0 


l'T„{y-)T^(y‘)dy 


»r, Nl Ki Ni Nj, . 

= -isi E E E E Ai4^N,^-^r7rfr^T^{y‘j)TM: 

■/V’ j=0fc=0n=0m=0 


hence 


ori ^y J' (y\)T fui) 

la. - l^ A:<PNxk c.a.c.G^ 

’ y j=0 fc=0 n=0 m=0 


where the integration term 6„m is defined in equation (2.46) 


(2.61) 
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In order to maximize /*, we require that the variations of this with respect to the 
nodal values on x = boundary vanish. In this maximization, 7’^ contributes 
only to those nodal equations which are obtained by differentiating w.r.t. nodal 
values which lie on the boundary x = a^. The contribution of 7’^ is 


„ N' Mi N< 

9 ^^y 

-jjT.Y.'LA. 

i=0n=0m=0 


TnmMv)) 


which is added to the right hand side of equation (2.38) for node 
Similar treatment can be done for other boundaries. 


(2.62) 


2.8 Spatial Derivatives at Collocation points 


Spectral element method allows extremely accurate evaluation of derivatives in each 
element in the computational domain. 

The spatial derivatives of the variable <f> with respect to x and y are given by 
equations (2.35) and (2.36). In element, equation (2.36) can be rewritten as 


d<t> 

dx 


9 9 Nx Ny jVi 1 

k=0 n=0 


Uii) 


dT^x) 

dx 


h{y) 


(2.63) 


Using the definition of Gauss-Lobatto Chebyshev points and that of Chebyshev poly- 
nomials, we have 


and 


Tnixj) = cos 



= COS 


riTtj 


dTnjx) 

dx 


n.sin [n.cos ^(x] 

\/(l 


(2.64) 


(2.65) 


substitution of (2.64) and (2.65) in (2.63) yields 


dx 


n n Nx Nx 1 

f|-EEE?>.w 

■Ltx -^^x A;=0 n=0 ^3 


rnrj ^n.sin[n.cos ^(x] 

— ^k{y) (2.66) 


N, 


VCl - 


As the function is represented by its values at the collocation points, its derivative 
will also be given by at these collocation points. For any point x = x^, we obtain, 
from equation (2.66) 


^ _ 2 2 

dXr 


^ Nx Ny Nx 1 

YJf S 

dyx jx^0k=0n=0 


.G. 


-cos( 


riTcj 


AT 




hk(y) 


(2.67) 
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For fixed y — yp, the term hk{y) in equation (2.67) is = 6kp and the summation over 
k can be dropped upon substituting <pjp for <f>jk. The resulting expression takes the 
form 


where 


and 


d(f) 

dXr 




yEKFniDP, 


jzzO n=0 


( 2 . 68 ) 


\ 

n.sm^ 

DPnr = - r^O or (2.70) 

sinj^ 

— r = 0 

= (-l)~-'n2 r = 


Similar treatment works to evaluate derivative with respect to y. 


Chapter 3 


Test Problems 


3.1 Unsteady Natural Convection in a Rectangu- 
lar Cavity 

The study of natural convection in enclosures is of considerable interest because of 
its importance in numerous engineering applications such as the cooling of electronic 
equipments, materials processing, solar energy systems, thermal energy storage sys- 
tems, and underground electric power transmission. Natural convection in enclosures 
is usually characterized by the formation of convective cells due to the buoyancy ef- 
fect. The heat transfer in enclosures is represented by correlating the average Nusselt 
number as a function of the Rayleigh number, Prandtl number, and aspect ratio. 
Enclosures can be classified in two categories: cavity-type and annulus-type. 

Natural convection in rectangular and square cavities has been extensively studied 
by a number of researchers [14- 18]. The problem has been tackled experimentally as 
well as numerically. Both finite-difference[15,16] and finite-element [17] techniques 
have been used for numerical simulations. 

To examine the usefulness and limitations of the present study, a problem of the 
natural convection in a rectangular cavity has been chosen. 

3.1.1 Physical Problem and Mathematical Model 

We consider a closed square two-dimensional cavity containing a Newtonian fluid 
(air) initially at rest at temperature Too ^ shown schematically in Figure 3.1. At 
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time t = 0, the left- and right-hand end walls are instantaneously heated and cooled 
respectively to temperatures and T^. and thereafter maintained at these tempera- 
tures. The upper and lower boundaries are insulated. All surfaces are rigid non-slip 
boundaries. The Boussinesq approximation is assumed to be valid, which assumes 
that fluid properties are constant everywhere and the density variation brings about 
the buoyancy term in momentum equations. The flow is assumed to be laminar and 
two-dimensional. The subsequent motion is described by the usual equations. 


3.1.2 Governing Equations for Laminar Flow 


For an incompressible viscous flow of Newtonian fluid the governing continuity and 

momentum equations can be written as : 

continuity 


du ^ 

dx dy 


(3.1) 


X momentum (note that x is the vertical direction) 


du du du 


dt 


dx ' dy‘ 


dp d^u 


y momentum 



dv dv. 


dp jd^v d^v^ 
dy ^ dx^ dy^ 


(3.2) 


(3.3) 


Energy 


dT dT dT ,d^T d^T, 

■ft + “ & + “ “'a? + 


(3.4) 


where u and v are the vertical and horizontal velocity components in the x and y 
directions respectively, p is the pressure, T is the temperature, p is the density of the 
fluid, p is the dynamic viscosity of the fluid, and a is the thermal diffusivity of the 
fluid. The viscous dissipation effect is neglected in energy equation. It is to be noted 
that g acts in the x direction. 

Following the Boussinesq approximation, density is everywhere constant except 
in the buoyancy term. The density variation with temperature is stated by the 
equation 


= p[i + a(r - T„)) 


(3,5) 
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Therefore 




(3.6) 


where subscript ‘oo’ stands for the reference state. The parameter 0 is the coefficient 
of thermal expansion of the fluid and is given by 

^ = (3.7) 

where T is the absolute temperature. For ideal gases, 0 = l/T. After substitution of 
(3.6) into (3.2) and simplification the x-momentum equation takes the form 


du , du 

p-^ + p{u 


dt 


du 

ax + = 


dp ,d^u a\ 


r„) 


(3.8) 


Here the term g0{T — Too) represents the buoyancy term. The velocity and the 
temperature fields are coupled. 

The boundary conditions for the present problem are 


cT 

II 

II 

f:=o 

OX 

on 

n = u = 0, 

^ = 0 
ox 

on 

II 

II 

II 

on 

u = V = 0, 

II 

on 


X = 0, 0 < y < L (3-9) 

X = H, 0 < y < L 
y = 0, 0 < X < H 

y = L, 0<x<H 


3.1.3 Non-dimensional Parameters 

The problem may be restated in the nondimensional form. We introduce here the 
following non-dimensional variables 


V=^ X=f 

Z) T-Tc p — (Poo-pW 

^ ~ Th-Tc ^ P‘'^ 

where i/ is the kinematic viscosity, Pr is the Prandtl number, and Gr is the Grashof 
number. It may be mentioned that in this problem Th has been considered as equiv- 
alent to of a general formulation. Similarly, instead of Too, we use Tc- 


U = f 

^ L2 


y-i 

Pr = !^ 

a 
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In dimensionless form, the equations are expressed as 

dll dv 


du , ^du ^du 


dr 


dx 


dX'^ dY ® 


^ ~ Tc)L^ 

' L o I 


dY dx ^dx^ 


av dP ,d^v d^v, 

dr ^ dx'^ ^dY dY ^dX^ ^ dY^^ 


(3.10) 

(3.11) 

(3.12) 


d9 de yde _ 1 d^e d^e 

' + r^,r — D- i V"2 


dr ' dX ' ^ dY Pr^dX^ ' dY^ 

The boundary conditions for the problem in nondimensional form are given by 


(3.13) 


u = 

II 

o 

^=0 

dX 

on X =0, o<y’<i 

u = 

II 

1^ = 0 
dX 

on X = Ar, 0 < F < 1 

u = 

II 

o 

0 = 1 

on y = 0, 0<A:<Ar 

u = 

II 

0 = 0 

on y = 1, 0 < X < Ar 


(3.14) 


where Ar = H/L is the aspect ratio. 

The conservation of the mass, momentum, and energy governing the buoyancy 
driven fluid flow and heat transfer in the enclosure may be formulated in terms of 
dimensionless quantities of stream function, vorticity, and temperature, respectively. 
For two-dimensional problem, it is convenient to define a stream function ^ and 
vorticity u as 




(3.15) 


" dX dY 


(3.16) 


We can readily see that existence of (3.15) automatically satisfies continuity equation 
(3.10). If we substitute the dependent variable with stream function, we shall not 
be concerned with equation (3.10) any more. Invoking equation (3.15) into (3.16) we 
obtain Poisson equation 
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Table 3.1: Expressions for parameters in general transport equation 


Function 


A 

r 


Temperature 

e 

1 

1 

Pr 

0 

Vorticity 

U) 

1 

1 


Stream function 


0 

1 

u? 


Now we differentiate equation (3.11) with respect to Y and equation (3.12) with 
respect to X. If we subtract differentiated equation (3.12) from differentiated equation 
(3.11) and rearrange the resulting equation, we shall obtain 


du) ^ fjdoj I ^ I 


(3.18) 


This equation is the vorticity transport equation. 

Instead of solving Poisson equation for stream function, we solve the Cauchy 
Kowaleska equation, i.e., the Poisson equation is recast in a false transient form as a 
solution strategy. 


6V’ 

dr 


+ "b ^ 


(3.19) 


dX^ ' 6F2 

The equations (3.13), (3.18), and (3.19) are the governing differential equations for 
two dimensional natural convection problems with moderate buoyancy effects. The 
vorticity equation is coupled to the parabolic Cauchy Kowaleska equation through 
the nonlinear convection terms and the energy equation is coupled to the vorticity 
equation. These dimensionless equations for vorticity, stream function and tempera- 
ture can be represented by the following general transport equation in terms of the 
flow property 




^^2 + 572!+ ^ 


(3.20) 


where is a source term, T is a diffusion coefficient, and ^ is a constant. The 
parameters used in general transport equation (3.20) are defined in Table 3.1. 

The boundary conditions to be satisfied by stream function, vorticity, and tem- 
perature are shown in the Table 3.2. 

The advantage of stream function vorticity formulation is that the pressure gets 
eliminated from momentum equations and we solve only two equations, one for stream 
function and other for vorticity, instead of solving three equations, one continuity and 
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Table 3.2: Expressions for boundary conditions for natural convection in rectangular 


cavity 


Function 

Temperature 

Vorticity 


Top wall Bottom wall Left wall Right wall 




e = i 


0 = 0 


Stream function V’ t/'or||- = 0 rp or ^ = 0 tp or M = 0 tp or ^ = 0 


two momentum equations. However, the approach^ot free from computational diffi- 
culties. In solution of vorticity equation, the first difficulty arises from the boundary 
conditions, since the vorticity is unknown along boundaries. The method generally 
employed is that the stream function equation is used at the wall in order to generate 
Dirichlet-type conditions. The second difficulty that arises is due to the variation of 
these boundary conditions at each time step. 

At no-slip boundaries, vorticity is produced. The boundary conditions for u are 
obtained using the no-slip condition at the wall and expression (3.17). Let us consider 
X = 0 wall, where 

We use finite-difference technique to compute vorticity at the wall. Applying central 
difference at a point ‘b’ on the wall, we get 

V’ui+l d” 1 /Q 

= 

2AX = 

where the suffix ‘ub’ stands for upper boundary. After simplification we get 


(3.21) 


(3.22) 


(3.23) 


In a similar fashion, we get 


X = Ar, 


r = o, 


y = i, 


= 


2V’66-1 

■ AX" 

' Ar" 

2V’r6-l 

■ Ar" 


(3.24) 
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subscripts ‘bb’, ‘lb’, and ‘rb’ represents bottom, left, and right walls respectively. 

The rate of heat transfer along each wall is determined from the temperature field. 
The local Nusselt number along the. heated wall is 


^r , 


(3.25) 


The average value of the Nusselt number is evaluated using 

1 fAr 00 

(3.26) 

The derivative ^ is computed by using spectral element technique. We compute the 
integration in equation (3.26) by the trapezoidal rule. 

1 N^e^+l 

£ [AXi{Nui-\-.Nui+^)] 


= ^ E [AX.-(iVu, + JVu,+,)] (3.27) 

where Nui is the Nusselt number at the point and AXi = — Xf. + 1 

represents the total number of points in x-direction. The domain is assumed to consist 
of tx elements in this direction. 


3.1.4 Discretization and Solution Procedure 

The governing equations are unsteady, coupled, non-linear partial differential equa- 
tions, for which a semi-implicit scheme is used for time discretization. As already 
stated in chapter 2, the non-linear convection terms are treated explicitly while the 
diffusion term implicitly. The time discretized form of the governing equations for 
temperature, stream function and vorticity is given by the well known Helmholtz 
equation of the form 

^xx + <i>yy ” = !<(> (3.28) 

Where the parameter and the forcing function for temperature, stream function 
and vorticity transport equations are defined in Table 3.3 

The spectral element method is used for spatial discretization, in which the vari- 
able (j) is expressed in terms of Lagrangian interpolant through Gauss-Lobatto Cheby- 
shev collocation points, which are closely spaced near the elemental boundaries. This 
method uses variational approach which requires that the functional is to be max- 
imized. Maxinoization of functional results in the elemental equations, defined in 
chapter 2. These equations are combined to form the global set which is expressed as 

/^global global _ global (3.29) 

where pg and Qg are the global indices, expressed by 
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Table 3.3: Expressions for and forcing function 


Equation (f> A^ 


Energy 

e 

El 

At 


Vorticity transport 

D 

1 

At 


Stream function 

i' 

1 

At 



Pg — 1)2,..., 4- l){Nyey + 1) 

9g — 1)2,..., {NxCx + l){Nyey + 1) 


and fq is the known forcing function which changes at each time step. The solution 
procedure used to obtain the steady state solution is described as follows: 

1. Specify the initial guess for the temperature, vorticity, stream function, and 
velocities. 

2. Solve the energy equation at n+1 time step. 

3. Solve the vorticity equation at n+1 time step. 

4. Solve the Cauchy Kowaleska equation for the stream function at n+1 time step. 

5. Compute the velocity field in the solution domain and vorticity at the bound- 
aries. 

6. Check the convergence by the criterion 


rms 


i 


- 'lif 

Ntotal 


< e 


(3.30) 


Where e is the preassigned tolerance value, and rms is the rms difference value, 
defined as the square root of the average of the squares of the differences between two 
successive time step values. 

The procedure is repeated until the convergence condition Eq.(3.30) is satisfied. 

The spectral element method has two ways to achieve better numerical accuracy, 
i.e., by increasing the nodes in the elements, and/or increasing the number of elements. 
The optimum choice of these two parameters depends on each individual problem to 
be solved. Without changing the total number of the interpolation points, using more 
elements in each direction reduces the computational cost, but at the same time, the 
accuracy level becomes lower. 
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3.2 Natural Convection Heat Transfer in a Verti- 
cal Rectangular Enclosure with Discrete Heat- 
ing 

The ever-increasing demand for faster and smaller computer elements has led to the 
development of compact electronic equipments of high heat dissipation rates per unit 
area [21]. Typically, the power density of ceramic chip carriers is about 500 KW/m^, 
and the maximum temperature allowed is approximately 100° C or less. For safe 
and reliable working of these densely packed electronic components proper cooling 
is essential to limit the operative temperature. Out of various cooling methods, 
natural cooling by free convection plays an important role in many electronic packages. 
Cooling by natural convection is preferred because of its high reliability, absence of 
noise, and low maintenance cost. Under certain circumstances, such as operation of 
electronic equipment in dusty or hazardous environment, electronic equipments are 
packed within sealed enclosures. The components may be mounted to one vertical 
wall of the enclosure, while one or more of the other walls is cooled. 

3.2.1 Physical Problem and Mathematical Model 

For the purpose of numerical study, we consider the natural convection cooling of 
discrete fluish-mounted heaters inside a vertical rectangular enclosure [19], as shown 
schematically in Figure 3.2. A vertical rectangular enclosure with four discrete fluish- 
mounted heat sources of isoflux qn on the left vertical wall is considered. The discrete 
heaters of length I are placed at an equal spacing s between the fixed adiabatic top 
and bottom surfaces of the enclosure. The right vertical wall of the enclosure is 
isothermally cooled at temperature Tc- 
The following assumptions are made: 

• Two-dimensional laminar flow and heat transfer. 

• The fluid (air) within the enclosure is considered to be Newtonian. 

• The physical properties, such as viscosity, thermal conductivity, etc., are as- 
sumed to be independent of temperature except for the density, in the buoy- 
ancy term in momentum equations, that is assumed temperature dependent. 
(Boussinesq approximation) 
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Figure 3.2: Physical configuration and coordinate system 
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With the assumptions stated above, the governing equations for buoyancy driven 
fluid flow and heat transfer in the enclosure, in terms of stream function, vorticity, 
and temperature are as follows 


dr dx dY Pr^dX^ 


(3.31) 


dr dX dY dY^^^ dY 

dxl^ d'^xj) di^xj} 

~ + 7:77:: + w 


dr dX^ dY^ 
with the boundary conditions: 

de 


t/> = 0, 

dY 


dY 




{at heaters) K = 0, 


= 0 


{elsexvhere) 


«- 

II 

0 

e=o Y=i 


de 

0 

11 



96 „ 

II 

0 



(3.32) 

(3.33) 

(3.34) 


where Ar is the aspect ratio (= H/L) is the aspect ratio, Gr * is the modified Grashof 
number defined as Gr* — 

The local heat transfer rates at the discretely heated wall and the isothermally 
cold wall are estimated by the knowledge of local Nusselt numbers defined as follows 


Nuh = ^ (3.35) 

Nuc = -(^)y=i (3.36) 
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3.3 Flow Analysis in Vertical Fin Arrays 

Vertical fins are widely used in the cooling of electronic devices. A numerical inves- 
tigation for predicting the heat transfer from the vertical fins have been reported by 
Sukhatme and Saikhedkar [22]. Present study uses a similar geometrical configua- 
ration for the Spectral Element Method. Essentially, the problem is concerned with 
vertical fin arrays which are attached on a hot base with temperature T^. The fins 
transfer heat to the surrounding fluid by natural convection. The vertical fin arrays 
to be analysed are shown in Figure 3.3. It consists of a large number of rectangular 
fins of a small thickness projected out from a common vertical base plate. Two ad- 
jacent fins and the base plate constitute a channel through which natural convective 
flow takes place. Due to the geometrical symmetry, the rectangular channel formed 
by the ABCD and EFGH planes form the domain of primary interest. However, 
because of the fact that the fins are situated in an otherwise quiescent environment 
which is influenced by the fin temperature, the analytical domain is extended and a 
rectangular channel between the planes KLMN and PQRS becomes the configuration 
of interest. The distribution of elements on each vertical plane is same and one such 
typical distribution on KLMN plane has been shown in Figure 3.4. The shaded area 
indicates the location of the fin. 


3.3.1 Physical Problem and Mathematical Model 

As mentioned earlier, a series of fins of height E in x-direction, length L in y-direction 
and the spacing 2B in z-direction are attached to the hot base which is maintained 
at uniform temperature Tu,. The solid surface and the fins constitute a part of the 
computational domain. Heat is transferred from the solid base of uniform temperature 
to fins by conduction and to infinite ambient fluid of temperature Too by natural 
convection. 

Let X denotes vertical coordinate direction, y denotes the coordinate perpendicular 
to the wall, and z coordinate is in the direction of successive fins. The gravitational 
field acts in x direction. The fluid flow is assumed to be essentially two dimensional 
with velocity components u and v being in x and y directions respectively. Cross flow 
in z direction is negligible in magnitude. Instead, it is assumed that the diffusion 
occurs in all the three directions. This simplified assumption of two dimensional 
flow allows us to use stream function vorticity formulation. Flow is considered to 
be laminar, further it is assumed that there is no viscous dissipation. The usual 




system. 
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Boussincscj approximation is assumed to be valid. The governing equations that 
describe this problem are the Navier-Stokes and energy equations. 

For an incompressible flow of a Newtonian fluid, the governing equations, incor- 
porating the above mentioned assumptions, are as follows: 
continuity 


du _ rv 

dx dy 


(3.37) 


X momentum 


du . du 



dp .dPu d^u 


di^u 

I? 


] - P 9 ^{T - 


Too) 


(3.38) 


y momentum 


dv dv dv. dp jd^v d^v d^v^ 

'“'a? + V + 


(3.39) 


Energy 


dT dT dT ,d^T d^T d'^T, 


(3.40) 


where u and v are the vertical and horizontal velocity components, Tko is the temper- 
ature of the quiescent ambient fluid, p is the dynamic viscosity and a is the thermal 
diffusivity of the fluid considered. Except for p in the buoyancy term in x momentum 
equation, the fluid properties are assumed constants. The dependence of density on 
temperature is related by 


^'oo = c[l + a{I’-r„)l (3.41) 

In the above equation, constant /? stands for the volume coefficient of thermal expan- 
sion of the fluid. 

3.3.2 Non-dimensional Parameters 

To represent the problem in non-dimensional form the following nondimensional pa- 
rameters are introduced 


U = ^ 


V=^ 


^ = f 


Y z= ^ 
^ H 


r = 


i/t 

H2 




Tui— Toe 


P = 


iPoo-p)H^ 


Z = W 




Pr=^ 
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Gr - 

i/2 


Here H is the lieiglit of the fin, is the reference pressure at reference temperature 
Too, i' is the kinematic viscosity and Pr is the Prandtl number of the fluid. Gth is 
the Grashof number based on the height of the fin. 

The governing equations for conservation of mass, momentum, and energy are 
nondimensionalized using these non-dimensional parameters. In non-dimensional 
form these equations are given by 


dX'^ dY~^ 


( 3 . 42 ) 


dr^ dx dY dX ^dx^ ^ 


d^U d^u 

QY 2 + 


-Grjjd 


( 3 . 43 ) 


dV dV dV dP d'^V d^V d'^V 
dr dX dY ~ dY'^ ^dX^ ^ dY^ dZ^^ 


( 3 . 44 ) 


do de do 1 dH dH d‘^e 

dr ^ dX ^ ^ dY ~ Pr^dX'^ ^ dY^ dZ^^ 


( 3 . 45 ) 


There are two ba.sic approaches for numerical solution of these equations: primitive 
variables approach and stream function vorticity approach. The primitive variables 
approach considers the above equations directly, taking the velocity components and 
the pressure. The stream function and vorticity approach employs derived variables, 
namely, stream function and vorticity to solve the problem. As a consequence the 
pressure gets eliminated and the existence of stream function identically satisfies the 
continuity equation. 

We have already defined vorticity w as 


- dX dY 


( 3 . 46 ) 


If we define stream function in two-dimensional plane as ^ — tp{X^Y), then 


U = 


dxjj 

dY 



( 3 . 47 ) 


which culminates in 


( 3 . 48 ) 
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The false transient form of equation (3.48) is given by 


di< dS> 

dr - ^ (3.49) 

Elimination of pressure from momentum equations results in the following vorticity 
transport equation 




(3.50) 


We note from the above equations that the Grashof and Prandtl numbers are the two 
governing parameters for the problem of interest. 

All these equations for vorticity, stream function, and temperature can be repre- 
sented by the general equation of the form 




d^<j> 


+ 


+ 


dX^ dY^ dZ^ 


r] + 


(3.51) 


where meaning of terms is explained in Table 3.1. 


3.3.3 Boundary Conditions for Present Problem 

The relevant boundary conditions for the problem are as follows: 


• At the fin base surface AEHD {2H < x < 3if,0 < 2 < = 0) (solid surface) 

The no slip boundary condition and the prescribed wall temperature is 

u = v = o i’ = o T = t ; 

dv A , , du d^rl / 

^ dy air 

• At the top extension of fin base surface AEPK (0 < x < 2H,0 ^ z B,y = 0) 
a; = 0 V- = 0 f = 0 

• At the bottom extension of fin base surface DHSN (3/f < x < AH, 0 < z < 
B,y = 0) 


• At the surface ABCD, z = 0 
symmetry conditions give 
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• At the vertical mid fin channel section PQRS 2 = 5 
symmetry conditions requires 


du A fh; rv 

dz ““ dz 

— = 0 

dz 

o 

II 

At the surface LQRM 



^ = 0 ^ = 0 

dy dy 

II 

o 


At the surface KLMNDCBA (z = 0) 


du A dv A 

dz dz ^ 

^ - 0 
dz 

— = 0 
dz 


The mass, momentum and energy conservation equations discussed above are the 
fundamental equations for natural convection in the geometry shown in Figure 3.3. 
The basic mechanism by which heat is transferred to the fins attached to surface 
is conduction. The fins are assumed to be very thin, equal to the thickness of the 
plane, and that the heat transfer from the fin boundaries is neglected. Instead, heat 
is transferred from the fin surfaces to ambient fluid in z direction. The temperature 
distribution in fin (solid) is governed by the transient heat conduction equation given 

by 


dZ _ kdT 

PsCs - K[ sdz 

where T, is the temperature in the fin (solid), /?,, c* and kg refer to the density, specific 
heat and thermal conductivity of the fin, S is the half thickness of the fin, and k is 
the thermal conductivity of the fluid. The second term on the right side of equation 
(3.52) expresses the heat loss from fin surface to the adjacent fluid on its either side. 
The boundary conditions along the fin boundaries can be taken as 


Tg = 25<x<35, y = 0 (3.53) 

BT 

^ = 0 on x = 25, 0<y<L 

ox 

^ = 0 on 25<x<35, 3 / = T 

ay 

^ = 0 on x = 35, 0<j/<L 

OX 


alongwith the initial condition Ts = Tso for < < 0. 

Since the energy equation (3.52) is essentially for the fin material, it is nondi- 
mensionalized separately. Here x, y, z and t are nondimensionalized in the same 
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way described earlier. Employing the nondimensionalization of T, in terms of refer- 
ence temperature 7U as 0. = (7'. - 7'^)/(7’, - r„), the governing equation (3.52) is 
obtained in diinensioiilesss terni.s, given by 

X pc Hde 

dr a Fr^dX^ Pr p,c, 6 ^ 

together with the initial and boundary conditions 


0. 

, = 1 

Os = O^o 

2 < < 3, 

for T <0 

y = 0 r > 0 

dO, 

Tx 

= 0 

on 

X = 2, 

0<Y<LIH 

dO, . 

Ih’ 

- 0 

on 

2 < a: < 3, y = L/if 

do. 

dX 

= 0 

on 

a: = 3, 

0<Y<L/H 


where a, and <v are the thermal difFusivities of solid and fluid respectively, v is the 
kinematic, viscosity and Pr is the Prandtl number of the fluid. dOJdZ is the temper- 
ature gradient in z direction. 

The above formulation suggests that the heated surface is cooled by free convection 
to the ambient fluid and by conducting heat to the fins which, in turn, convect heat 
to the ambient fluid. The convective heat transfer coeflBcient can be presented in 
terms of the local Nusselt number Nu, which is defined as Nu = Also the local 
heat transfer coefficient h may be expressed in terms of the temperature gradient at 
the surface, ambient temperature and surface temperature. Thus 


Nu = 


k 



(3.56) 


Here, the subscript w represents the relevant wall or surface, k stands for the thermal 
conductivity of the fluid and is the temperature gradient normal to the surface. 
The mean Nusselt number Nu for the surface is defined as Nu = where 


Nu=\ [ NudA = i / Nuds (3.57) 

A JA o Js 

In above relation for average Nusselt number, ‘A’ is the total surface area exposed to 
the fluid and S is the distance along the surface. 

The local Nusselt number at the fin surface, in terms of non dimensional quantities, 
is expressed as 


de 

Nu{X,Y)^-^ 


(3.58) 
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We us(- the serond-orcler Taylor expansion to approximate the derivative of 6 on the 
surface, d’he resulting expression is given by 


;vu(.Y,r) = 

2AZ 


(3.59) 


The spanwise average Nusselt number distribution on the surface, Nu,^x^ is given by 


] rL/H 

NviS,X) = — Nu{X,Y)dY 


(3.60) 


Ny.ey+l 


LjH 

~ 2Lj H ^ ■i' 

and the area averaged Nusselt number, based on the area of the fin surface, is obtained 
from 


yVu(S,a) = C Nu{s,X)dX 
Jo 


(3.61) 


1 iVx.ei-f-l 

- £ [AJi:i(7Ki+v;,-„)] 

^ t=l 


3.3.4 Discretization and Solution Procedure 

Solution for a three-dimensional flow problem requires considerable effort. However, 
with some simplifications the problem can be dealt easily. 

Comparison with the two-dimensional problem shows the presence of one addi- 
tional diffusion term in governing equations. This term, along with non linear 
convective and source terms, is treated explicitly and the two-dimensional diffusion 
in X and y directions implicitly. This treatment enables us to adopt the same spatial 
discretization discussed in chapter 2 with only one difference of inclusion of this addi- 
tional term in forcing function. To simplify the analysis further, the finite difference 
scheme is used to approximate all the derivatives in z direction and the problem is 
solved for successive planes in this direction. 

The time discretized form of the equation (3.51) is given by 


4' 




4^ A 


+ v?± 

dX^ dY 


d^4 av 
dx^ ay2 


n-t-l 


-h 




dZ^ 


+ 

qn 


(3.62) 


this can be rearranged as 

Vax^ ) + ) 




<t>^ A 

tat r 


dx dY 




SI 


(3.63) 
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'ral»!(‘ 3.4: Exi)r(>ssi()nK for and forcing function 


K<juat ion 

<i> 


U 

Energy ( 1 ) 

0. 

(Kpr 

At 

otPrd'^ 1 ct pc H / dd yi 
q$At ' Os PsCs 6 

Energy (2) 

0 

El 

At 

+ (Sr 

Vorticity transjuirt 

u.’ 

1 

At 

-si + Pli + vf 1" - Gr(f r - (&)" 

Stream function 

i' 

1 

At 



This is th(' Hclrnhoh// <‘quation of the form 

<t>Tx + 4>yy ~ — fd) (3.64) 

where Aj = 1/(1 ' At) is a parameter that can take different positive values and ^ 
is known forcing function that forms the right side of equation (3.64). Equation 
(3.54) can be discredited in the same way, treating the diffusion term implicitly and 
source term exjilicitly, and tlie discretized form can be also be represented by equation 
(3.64). Definitions of the various terms that describe the energy, stream functon, and 
vorticity equatiiins are summarized in Fable 3.4. 

The central diff<‘r(>ncing of the diffusion term yields 

d'^(t> _ <l>{k + \) - 2(l){k) + <t>{k - 1) 

and the first ordt'r derivative in source term of equation (3.54) can be differentiated 
using forward diff<‘renc.(' to obtain the following approximation: 

^ _ ^(fe + 1) - <f>{k) (3.66) 

dZ ' ' 

Numerical solutions to th<‘ foregoing governing differential equations for the present 
problem ar(* obtained by spectral element method. The spectral element discretization 
scheme used in this problem is presented in chapter 2, and the solution methodology 
employed here is basically the same as for rectangular enclosure problem with one 
extra step of marching in z direction after having obtained the steady state solu ion 

for previous z-plane. 



Chapter 4 


Results and Discussion 


The computat ions have l)e<>n accomplislied for buoyancy driven j9ows and heat transfer 
in rectangular tuiclostires and in vertical fin arrays. Results, so obtained, have been 
compared witli the finite-difF<‘r<-nce and finite-element solutions to assess the accuracy 
of th(‘ spectral ehnnent method. 

Tlie num<‘rical work has I teen performed on the HP 9000/735 series computer. 
The system has a RAM of 144 MB and 40 Mflops (scalar) peak performance. 


4.1 Unsteady Natural Convection in a Square Cav- 
ity 

Calculations have heim performed considering air as the working fluid (Pr=0.71) for 
aspect ratio Ar=LU and 10'* < Ra < 10“*. The flow and temperature fields inside 
the cavity are illustratt'd by rmmns of contour plots of streamlines, isotherms and 
vorticity, respectively, as shown in Figure 4.1 to 4.9. The results shown in these 
figures are t)l>tained for dinx'nsionless left side wall temperature equal to 1.0 and 
that of right side wall temperature equal to 0.0. The top and bottom walls are held 
insulated. Th<* initial tiunperature, vorticity and stream function fields are chosen 
to be zero at. time t = 0. Results are obtained using 12 x 12 Gauss-Lobatto points 
and one element in each direction. Similarly 6x6 and 4x4 Gauss-Lobatto points 
for each ehunents are used while the number of elements in each direction are 2 and 
3 respectively. It was observed tliat a system of one element in each direction with 
12x12 Gauss-Lobatto points is a good choice. The agreement between the present 
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Table 4.1: Average Nusselt number along the heated wall for air as working fluid (Pr 
= 0.71) Ar= 1.0. 


Ra 

Present study 

Philips [15] 

CO 

O O 

1.16438 

2.24574 

1.118 

2.250 


Table 4.2: Avera ge Nusselt number along the heated wall for Pr = 1.0, Ar = 1.0 . 


Ra 

Present study 

Shiralkar and Tien [14] 

10^ 

1.068 

1.131 

lO'* 

2.2774 

2.277 


spectral element solution and the finite difference solutions by Philips [15], Kublbeck 
et al. [16] and finite element solution by Ramaswamy [17] is excellent in all respects. 
The comparisons have been done on the basis of streamlines, isotherms and vorticity 
contours. 

For small values of Rayleigh number (< 10^), the isotherms are evenly spaced 
parallel lines indicating that the heat transfer is dominated by conduction and the 
convection effects are less (see Figure 4.1). For higher Rayleigh numbers (> 10 ), 
convection effect starts dominating as it is evident from distorted isotherm patterns 
shown in Figure 4.4 and 4.7. Isotherms for higher Rayleigh numbers result in high 
temperature gradients near the wall. 

Figures 4.2, 4.3, 4.5, 4.6, 4.8 and 4.9 show streamlines and constant vorticity 
contours in the cavity for three different Rayleigh numbers. 

The local Nusselt number distribution along the heated wall of the cavity is re- 
ported in Figure 4.10 for Ra = 10^ and 10^. The Nusselt number gradually decreases 
along the vertical wall. It is interesting to note from this figure that the maximum 
value of the Nusselt number occurs at a location slightly above the lower wall. 

The average Nusselt number along the heated wall is documented in Table 4.1. 
Results are compared with those of Philips [15]. Philips [15] had reported the results 
due to finite difference method using 65 x 65 grid pints. The results for the Spec- 
tral Element method using 12 x 12 Gauss-Lobatto points show a reasonably good 
agreement with the reported results. 

Again for the purpose of comparison of our results with those of Shiralkar and 
Tien [14], the same problem has been solved for Pr=1.0 and for a range of Rayleigh 
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number. The Nusselt number values due to present computation are reported in Table 
4.2. Table 4.2 also shows the corresponding Nusselt numbers predicted by Shiralkar 
and Tien. The results of present study agree well with the results of Shiralkar and 
Tien. 

It has been observed that when Rayleigh number increases and becomes equal or 
greater than 10®, the solution starts oscillating and the steady state values are not 
obtainable. It is conjectured that the flow displays an oscillating behaviour which 
corresponds to a typical time-periodic motion. 
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4.1.1 The Effect of Aspect Ratio 

To study th(' (‘frcct of a,si)ect ratio on the temperature and flow field, a problem of 
natural (■oiiv<Ttion in n'ctangular enclosure is solved. The temperature boundary 
conditions taken on tlx’ left and right side walls are 0.5 and -0.5 respectively with 
adiabatic top and bottom walls. Other boundary conditions remaining the same as 
explained earlier. Th<' initial conditions for 0, w, and V’ at r = 0 are zero. The 
calculations hav(‘ Ix't'ii performed for Ra = 1.446 x 10'* and Pr = 0.733. Figure 4.11 
shows the isotlu'rms, streamlines and vorticity contours for a rectangular enclosure 
of aspect ratio 3. Idie present results agree qualitatively with those of Ramaswamy 
117], 

Finally, Figures 4.12 and 4.13 show similar results for a rectangular enclosure of 
aspect ratio 1 .83, with the temperature boundary conditions -0.5 and 0.5 respectively 
on left and right side walls. The results shown in these figures are obtained for a 
Rayleigh number of 10^ and Pr = 1.0. Results for the same problem of aspect ratio 
of 1.83 but with linear temperature distribution on the horizontal walls are reported 
in Figures 4.14 and 4.15 . The top boundary temperature distribution is governed by 
the linear relation 0(0, 7) = -0.5 + Y while the bottom boundary of the enclosure is 
maintained at 0(1, F) = -0.4. Szekely and Todd [19] have reported the experimental 
and finite difference solution results for this problem for Ra = 8200 and Pr = 
2450. 

It was observed that for larger aspect ratios the numerical solutions converged 


slower. 
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4.2 Natural Convection Heat Transfer in a Verti- 
cal Rectangular Enclosure with Discrete Heat- 
ing 

A coiTiputational study of natural convection in an enclosure as applied to applications 
in cooling of electronic components is reported. 

The computations are first carried out for aspect ratio (Ar = 6), Ra* = 10® and 
Pr = 0.71. The relative heater size, i/H, and the location, s/H, on the vertical wall 
are fixed at 1/30 and 13/75, respectively. The qualitative behaviour of streamlines 
and isotherms inside the discretely heated enclosure is found to have good agreement 
with the results of Ho and Chang [20]. Ho and Chang used a mesh system consisting 
of a uniform grid of 151 in the x direction and a non uniform grid of 51 in y direction 
to generate the numerical simulations for the present problem. The present analysis 
employs 9 elements in x-direction and 1 element in y-direction. In each element 7 
Gauss-Lobatto points have been used in x- and y-directions respectively. 

At Ra* = 10®, heat of the discrete heaters is essentially dissipated via a conduction- 
dominated mechanism as indicated by the isotherm pattern around the discrete 
heaters shown in Figure 4.16. With increase of the aspect ratio of the enclosure, 
the buoyant convection flow is increasingly strengthened, exhibiting a transformation 
from a primarily clockwise unicellular flow into a structure characterized by regularly 
spaced multicellular flow in the core region for Ar > 6, as depictedby the streamlines 
in Figure 4.17. These split core vortices are essentially in parallel with the location 
of the heaters on the vertical wall. 

Figure 4.18 shows the predicted isotherms and streamlines pattern for i?o* = 
4.0 X 10^ and Ar =10. It is noteworthy that the present predictions agree qualitatively 
well with the experimental and numerical predictions of Ho and Chang [20]. 
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Figure 4.16: (a) Isotherms, (b) streamlines and (c) vorticity contours for discretely 
heated rectangular enclosure at Ra* = 10^, Pr=0.71 and Ar=6. 
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4.3 Flow Analysis in Vertical Fin Arrays 

I Ilf foinimtiitioiis have lircii carried out assuming air (Pr=0.71) as the working fluid. 
The results are reported for Rayleigh numbers of 10'* and 10®. 

The temperature distribution on the fin surface is obtained by solving the transient 
heat ronduetion etiuation. Subsequently the advection diffusion equation is dealt 
with. The fins are assumed to have a high thermal conductivity. The following 
property values have been used: 

O', = 8.418 X 10-® mV.s, p, = 2707 Kg/m^, jk, = 204 W/m°C, 

Cp, = 0.896 KJ/Kg°C 

As mentioned earlier, air has been assumed as the working fluid and the following 
property values have been used for air wherever needed 

= 0.2216 X 10"“ rnV^, Fa = 15.69 x 10"® mV^, Jka = 0.02624 W/m “C, 

The fins are jutting out from a vertical base. As they have negligible thickness; 
the H/S is as high as 10 in the present computations. The Rayleigh numbers used 
in the present problem are 10“ and 10®. The Spectral Element technique has been 
deployed to solve the governing equations for vorticity transport. Fig.3.4 shows 4x4 
elements on the x-y plane. The fin is an element which has a count as the third 
element in x-direction and the first element in y-direction. The height and the length 
of the fin are taken as 1 .0 ( nondimensional ). 

Figures 4.19(a) to 4.25(a) show the isotherms on different x-y planes at an interval 
of AZ = 0.16 for a Rayleigh number (Ra) of 10“ and Prandtl number (Pr) of 0.71. The 
first x-y plane [Fig. 4.19(a)] shows the isotherms on the fin element and the immediate 
neighbourhood. Figures 4.20(a) to 4.25(a) depict the temperature distribution of 
the flowing fluid at different vertical planes. The temperature distribution at the 
symmetric mid-plane is shown in Fig. 4.25(a). It can be seen that the temperature 
is advected markedly in the x-direction. 

Figures 4.26(a) to 4.32(a) show the isotherms on various x-y planes for a Rayleigh 
number of 10®. The same qualitative trend of distribution of the isotherms, as it has 
been seen for Ra = 10“, is observed in all the cases. However, the quantitative values 
of the isotherms for Ra = 10® are somewhat different than those for Ra = 10“ on the 
same x-y planes. 

Figures 4.19(b) to 4.25(b) show the streamlines on seven different x-y planes at an 
interval of AZ = 0.16 for a Rayleigh number of 10“. A buoyancy driven flow pattern 
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'I'ablc' 4.3; Area a verage Nusselt numb er on fin surface. 
_Ra 10'* 10® 

Pr 0.71 0.71 
^.,a 8.99 9.23 


is discerned IhroTigh the numerical flow visualization on each x-y plane [essentially 
Fig. 4.20(b) to 4.25(b)] mentioned above. The streamlines are formed by the cold 
air entering from the bottom and front openings through multiple plumes. This flow 
structure basically follows a single-chimney like flow pattern. The flow structure 
is self-sustained since fin to air heat transfer promotes the fluid motion. Figures 
4.26(b) to 4.32(b) illustrate the streamlines on all seven x-y planes, eis described 
earlier, for a Rayleigh number of 10®. The streamline distribution on all the planes 
confirms a single- clfmney type flow structure through the arrays. It is evident that 
with increasing Rayleigh number (so to speak, with increasing temperature difference) 
the cold air enters from the open end at the bottom and leaves through the open end 
at the top. 

Due to pausity of time, a detailed parametric study for various Rayleigh number 
and Prandtl number could not be accomplished. However, for two different Rayleigh 
numbers of 10^ and 10® the values of area average Nusselt number are reported in 
Table 4.3. 

For a similar geometrical configuration like the present investigation and for a 
Rayleigh number of 2.593 x 10“*, Shalaby et al. [23] obtained a value for Nu^a as 
9.05. This does not directly indicate the accuracy of our computation but confirms 
the acceptability of the predictive procedure. 

The value of Tfu^a for a Rayleigh number of 0.71 x 10® on a vertical isothermal 
plate of height 10 cm losing heat by natural convection comes out be approximately 
20. Thus it is seen that the values for the array are lower than for a plate alone. It 
will be interesting to investigate whether the Nusa value for the fin approaches the 
Nusa value of the plate as the spacing increases for an identical Rayleigh number. 
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Chapter 5 


Conclusions 


P„U«n. of buoy.nc. driven fluid flows and su— 

the problems is seim-implicit, the time stepp g 
obtained between the results of earlier ^ 

For the problem of buoyancy of Ho and Chang [20]. 

rr^ r :r:C:-;:a ar^t parses .r computers and 

other electronic applications. r i -irarirv driven flow in a vertical 

Finally sorne useful ^ 2: of the heat 

tilt ar: plcted. The results are compared with the previous numer- 
ical results based on simplifying assumptions. 
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